Studies of non-magnetic impurities in the spin-1/2 Kagome Antiferromagnet 
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Motivated by recent experiments on ZnCu3(OH)6Cl2, we study the inhomogeneous Knight shifts 
(local susceptibilities) of spin 1/2 Kagome antiferromagnet in the presence of nonmagnetic impuri- 
ties. Using high temperature series expansion, we calculate the local susceptibility and its histogram 
down to about T = 0.4J. At low temperatures, we explore a Dirac spin liquid proposal and calculate 
the local susceptibility in the mean field and beyond mean field using Gutzwiller projection, finding 
the overall picture to be consistent with the NMR experiments. 



I. INTRODUCTION 

Spin liquids are some of the most interesting spin 
phases that we know to exist theoretically. However 
experimentally they are hard to achieve because of the 
competition with magnetically ordered or spin-Pcicrls 
phases. Frustration can suppress the order, and one of 
the promising systems where strong frustration occurs 
with natural interactions is the nearest neighbor spin 1/2 
antiferromagnet on Kagome lattice. 

Herbertsmithitc ZnCu3(OH)6Cl2 contains Kagome 
layers of spin-1/2 magnetic Cu^"'' moments and 
has emerged as one such candidate in the recent 
years. i^'^'^'^'^'^i'^ Despite the large Curie- Weiss temper- 
ature of roughly 300K, experiments observe no magnetic 
order down to 50mK and no sign of spin freezing or spin- 
Peierls transition. Various measurements also suggest 
that there is no gap to excitations, but the nature of the 
phase remains elusive. 

From the beginning^i^ it has been appreciated 
that there is significant amount of disorder in the 
ZnCu3(OH)gCl2. Even though impurities do not seem 
to produce a glassy behavior, they can mask the intrinsic 
thermodynamic properties of the system such as the spin 
susceptibility and specific heat. The type of disorder has 
been clarified to some degree in recent experiments j^ii 
which suggest that 5 to 10% of the Cu that would ideally 
be in the Kagome layers interchange their position with 
the nonmagnetic Zn between the layers. The picture that 
is suggested is that the misplaced Cu behave as weakly 
coupled spins giving rise to the Curie tail and anomalous 
low-temperature specific heat in the absence of the field, 
while the misplaced Zn effectively create spin vacancies 
in the Kagome layers. 

There have been many theoretical 
studies of the clean spin-1/2 Kagome 
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Exact diagonalization studies^ i^^'^^i^^'^'^i"'^^ indicate no 
magnetic or spin-Peierls order. Spin wave and series 
expansion studies^iii^ near several ordered states also 
indicate that the magnetic orders are not stable (but a 
recent work2^ suggests that a valence bond solid may 
be stable). High temperature expansion studiesi^iii 
for susceptibility of the Heisenberg model provide good 
quantitative comparison with the experiments and 



recent worksi^ii^ account for possible Curie impurity 
spin contribution or possible Dzyaloshinsky-Moriya 
interactions to better match the measurements in the 
Herbertsmithitc. 

There arc a number of theoretical proposals for spin 
liquids in the spin-1/2 Kagome ^^^^a^iSS^^^ Of 
interest to the present work arc gaplcss spin liquids. 
One is the Algebraic (Dirac) Spin Liquid (ASL) of 
Hastings2i and Ran et ali^ whose trial energy is very 
close to the ground state energy obtained from the ex- 
act diagonalizatioiii^ Recently, Ma and Marston^i sug- 
gested that a different spin liquid with a spinon Fermi 
surface may be stabilized in the presence of ferromag- 
netic second-neighbor interactions. Coming from a very 
different perspective, Ryu et al^ suggested an Algebraic 
(Dirac) Vortex Liquid, which bears some resemblance to 
the ASL in its gaplcss nature but is a distinct phase. 

Impurities modify the properties of the ground state, 
but at the same time they can be used as a probe of the 
system, which is our pursuit here. In the NMR context, 
impurities allow the external magnetic field, which would 
be strictly q = probe in a clean system, to couple to 
other potentially more important wavevectors. Despite a 
large body of work on defects in correlated systems (see 
Refs. [30II31I and references therein), there have been only 
few studies of impurities in frustrated magnets ^ii^^i^i^i 
Exact diagonalization study by Dommange et al^ pro- 
vides important hints on the behavior of spins near a 
vacancy in the Kagome antiferromagnet. 

Our work is motivated by the recent ^^O NMR mea- 
surements of Olariu et al^, which shows direct access to 
the local spin susceptibility. Each oxygen atom feels the 
magnetic moments of two copper atoms it is connected 
to, see pictures in Refs. HQ- Because of the nonmagnetic 
impurities, some oxygens feel two and some only one spin. 
The two peaks that are observed in the NMR spectra at 
high temperatures are explained as coming from these 
two types of oxygens, neglecting variations in local sus- 
ceptibility. As the temperature is lowered the peaks first 
move towards larger susceptibilities and then back to the 
smaller ones broadening at the same time. At the lowest 
temperatures only one peak is visible which saturates at 
small but positive susceptibility value. A strong inhomo- 
geneous broadening was also observed in '^^Cl NMR at 
low temperatures^ but with no clear features that would 
allow spatial resolution. 
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In this paper we calculate local susceptibility of 
Kagome lattice in the presence of nonmagnetic impuri- 
ties, which are treated as missing spins, by two different 
methods. 

In the first part we use 12-th order high temperature 
series expansion for the nearest neighbor antiferromag- 
netic Hcisenberg model. We calculate local susceptibility 
near a single impurity Fig. [T] and for a lattice with 5% 
of impurities Figs. [2l [3] We obtain quantitatively ac- 
curate results down to about T k, 0.4J. These can be 
directly used to test whether the Heisenberg model with 
spin vacancies applies to the ZnCu3(OII)6Cl2. For this 
range of temperatures, the local susceptibility goes to the 
bulk value already at a few lattice spacings from an im- 
purity. On the other hand, the local susceptibility Xi 
of the spin next to the impurity is larger than the bulk 
Xbuik by as much as 15-20%. This can be understood by 
realizing that the Cu near the impurity has only three 
neighbors instead of four, which directly affects the local 
Curies- Weiss temperature giving 14% change in x this 
temperature. The series expansion study develops this 
observation systematically to significantly lower temper- 
atures and for all sites. On the other hand, we find that 
Xi starts to decrease around T k, 0.6J and drops below 
the bulk susceptibility at T w 0.4J, well before the bulk 
susceptibility starts to decrease at about J/6J^J^ This 
is consistent with the ground state picture found in the 
exact diagonalization study^^^ where a pair of spins that 
are both next to impurity forms an almost perfect singlet, 
so these spins become effectively non-responsive. 

One lesson we learn is that even at elevated temper- 
atures one cannot assume the same susceptibility on all 
Kagome sites. When connecting with the experiments, 
we suggest to consider at least two different susceptibili- 
ties: xi for the immediate neighbors of impurities and 
Xbuik for the rest of the Cu sites. In the ^^O NMR 
experiments)^ we then suggest four main groups of the 
oxygen Knight shifts as shown in Fig. [31 We hope that 
our finding and improved experimental resolution can be 
used to further test the applicability of the Heisenberg 
model to the Herbertsmithite. 

In the low temperature study, we attempt to un- 
derstand the susceptibility histograms in the spin liq- 
uid framework, taking up the ASL of Hastings^! and 
Ran et al^ with Dirac spectrum for the spinous. First 
we perform the mean field analysis treating the spinous 
as free fermions hopping on Kagome lattice with impurity 
sites removed. Again we calculate the susceptibility near 
a single impurity and in the presence of 5% of impuri- 
ties. We should say from the outset that such a simplified 
treatment of the spinous is likely inadequate near impuri- 
ties and one should allow modifying both spinon hopping 
amplitudes and fluxes near impurities. For example, the 
exact diagonalization study^^ finds strong singlet correla- 
tions on bonds next to an impurity. Furthermore, think- 
ing about possible more accurate self-consistent treat- 
ments suggests that the spinon hopping should also be 
strengthened on these bonds. Such modifications will 



certainly have significant effect on the calculated suscep- 
tibility of the spins nearest to the impurity. However we 
argue that there are robust features in the susceptibility 
histograms that are insensitive to the details of the treat- 
ment near the impurities. The vacancies are a source of 
disorder which will induce a finite density of states at 
the Fermi level with a characteristic spatial distribution. 
The bulk susceptibility associated with the sites that arc 
further from impurities decreases as the temperature is 
lowered, just as it would for the pure system where it 
goes all the way to zero because of the vanishing den- 
sity of states at the Fermi energy in the Dirac spectrum. 
However the impurities cause finite density of states and 
as more and more sites start to feel the impurities as we 
lower the temperature, which all sites eventually do since 
there is a finite density of impurities, the movement of 
the histogram peak to zero stops as is seen in Fig. [T] 
There, the detailed features seen at high temperatures 
depend on the treatment of impurities and should not be 
taken quantitatively; on the other hand, the overall pic- 
ture of the disorder-induced broadening of the vanishing 
pure Dirac susceptibility into an asymmetric histogram 
should be more robust. 

We also test the approach beyond the mean field by 
studying the magnetization distribution in the trial spin 
wavefunctions obtained by Gutzwiller projection. This 
approach can reproduce nontrivial results for the local 
susceptibilities near impurities in one dimensioui^ In the 
Kagome system with impurities, the main effect that we 
find is the increase in local susceptibility variation by 
about a factor of two, while the qualitative predictions 
of the mean field remain unchanged. In particular, we 
find that the local susceptibility can be also negative, 
deviating significantly from the positive bulk value. This 
can happen because there are strong antiferromagnetic 
correlations in the system. 

From the above discussion, it is clear that the local sus- 
ceptibility variation due to impurities is significant both 
at high and low temperatures and its experimental ob- 
servation can reveal a wealth of information about the 
spin state. At high temperatures, it is possible to track 
several features associated with different locations rela- 
tive to impurities. In particular, the susceptibility at the 
oxygen sites next to impurity should be xi which is dif- 
ferent from the half of 2xbuik seen by the oxygens in the 
bulk. If one could improve NMR line resolution, more 
peaks should be observable at high temperatures. On 
the other hand, the peaks broaden as the temperature is 
lowered and particularly as the system enters the corre- 
lated paramagnet regime. Our spin liquid study suggests 
that there remain no features in the susceptibility his- 
togram except for the overall peak associated with the 
bulk sites, which first moves towards zero as in the pure 
Dirac spin liquid and then stops because of the disorder- 
generated nonzero density of states at the Fermi level. 
This is reminiscent of what is observed in the "'^^O NMR 
experiments, even if not in every detail; it is also qualita- 
tively different from the spin liquid with a Fermi surface, 
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where there is a finite density of states in the pure sys- 
tem and the impurities would only broaden the histogram 
both ways without changing significantly the bulk peak 
location. 



II. HIGH TEMPERATURE SERIES 
EXPANSION 



temperature they start diverging and we take that as a 
point where the approximation stops being valid. Differ- 
ent values of a are tried, and sometimes it is possible to 
tune to a value where all the curves overlap completely to 
a much lower temperature, but that is a pathology, prob- 
ably indicating that the polynomials are all the same. For 
other values of a the curves usually start to diverge from 
each other at around the same temperature. 



We consider spin 1/2 nearest neighbor anti- 
ferromagnetic Heiscnberg model on Kagome lattice. 
Local spin susceptibility at site i is given by 



Xioc(i) 



{g^,Br{SfSU) 
knT 



(1) 



where S^^^ = J2j^j- We calculate x's in the high 
temperature series expansions in the presence of a sin- 
gle nonmagnetic impurity, treated as a missing site (va- 
cancy). We also study the Kagome system in the pres- 
ence of a finite density of nonmagnetic impurities, which 
we chose to be 5% as motivated by cxperimentsSiS'i^ in 
ZnCu3(OH)6Cl2. 

The expansion is performed to the 12-th order in J (or 
13-th order in 1/T) using the linked cluster expansion»^ 
The outline of the procedure is as follows. We generate 
all abstract graphs up to desired size. Then we gener- 
ate all subgraphs of these graphs, keeping track of the 
location of each subgraph in the graph. We calculate 
the local susceptibility for each graph at each point of 
the graph. Then we subtract all the subgraphs of each 
graph as needed in the linked cluster expansion to get the 
contribution this graph would give when embedded into 
lattice. The local susceptibility on any lattice can be cal- 
culated by creating all possible embeddings of all graphs 
and adding their contributions at every site. In this gen- 
eral formulation, the lattice does not need to be regular; 
in particular, this procedure applies to the Kagome lat- 
tice with impurities. In the practical implementation, we 
use the symmetries of the underlying pure Kagome lattice 
for efficiency and discard the embeddings that contain an 
impurity. In this way, we obtain exact 1 /T scries for the 
system with impurities. 

One more remark about the method. In the deriva- 
tion of the linked cluster expansion for the local suscep- 
tibility, one needs to calculate expressions of the form 
TT{H''SfH''S^^^H'') where a, 6, c are integers. Since the 
total spin is conserved, [S^^^,H] = 0, and because of the 
cyclic property of the trace, this can always be written as 
TT{H°'~^''~^'^Sf S^^^). Thus every diagram can be evaluated 
using this expression. 

After obtaining the series, we extend it beyond the 
radius of convergence using the method of Fade approxi- 
mants. We use [5,6], [5,7], [6,6], [6,7] and expand in vari- 
able l/{T + a) where a is usually 0.08 as in Ref . l37l. De- 
pending on a one might get a pole in the expression, and 
hence divergence in susceptibility even at relatively large 
temperature. This usually happens say in one of the ap- 
proximants while the others still overlap. At low enough 



A. Single Impurity 

The series coefficients of the pure Kagome lattice sus- 
ceptibility and of the local susceptibilities of the first 
four nearest neighbors near a single impurity are in Ta- 
ble [H The corresponding susceptibilities are plotted in 
Figure [T] 
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FIG. 1: [Color Online] Susceptibility of the pure Kagome sys- 
tem (far from any impurities) and local susceptibilities at the 
first four inequivalent neighbors of a nonmagnetic impurity. 
For each Xi, there are actually four curves of the same color 
corresponding to the Fade approximants [5,6], [5,7], [6,6], 
[6,7]; the plotting is cut where these curves start to diverge 
from each other. The X2 and xs happen to overlap over the 
entire temperature range shown, x is in units (qhb)^ /J- 



From Fig. [H we see that for the temperatures T > 
0.3 J, only the first few neighbors - most visibly the xi 
right next to the impurity - deviate significantly from the 
susceptibility of the clean system. The xi is about 15- 
20% larger than Xuniform over a wide temperature range 
O-^J S T S As the temperature is lowered further, 
the xi starts to decrease while the Xuniform still grows, 
and we see xi dropping below Xuniform- This is reminis- 
cent of the picture in the exact diagonalization study in 
Ref. m, in which a pair of spins in the same triangle as 
the impurity (there are two such triangles) forms a sin- 
glet at low temperatures, and therefore the susceptibility 
of these spins is strongly reduced. 
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TABLE I: Series coefRcients a„ of x = ^ X^n=o 4ra+i(n+i)! ( j')" foi' susceptibility Xuniform of the pure Kagome lattice and for 
the four closest neighbors Xi of nonmagnetic impurity as indicated in Figure [T] The coefficients for the uniform susceptibility 
agree with Ref. [l^ except for n = 6, which is probably a typo there. 



n 


Xuniform 


Xi 


X2 


X3 


Xi 





1 


1 


1 


1 


1 


1 


-8 


-6 


-8 


-8 


-8 


2 


48 


30 


42 


42 


48 


3 


-96 


-120 


-8 


-8 


-48 


4 


-320 


680 


-360 


-240 


-1360 





-38(84 


-14592 


A ccoc; A 


AT AO A 

-4(424 


-4U0i2 


6 


677824 


158144 


540624 


465248 


995680 


7 


14176256 


2843520 


16451456 


18525184 


15092736 


8 


-429202944 


-57813120 


-312132096 


-249716736 


-636643584 


9 


-11927946240 


-2620491520 


-12841824000 


-15158219520 


-12108702720 


10 


501186768896 


64507410176 


346403354880 


271966979328 


700664177664 


11 


13960931721216 


3019001330688 


14456058461184 


17966636077056 


13440259983360 


12 


-841802086780928 


-98970319958016 


-558124102202368 


-437482821907456 


-1111015315910656 



B. 5% Density of Impurities 



We take a 60 x 60 x 3 Kagome lattice with periodic 
boundary conditions and 540 randomly placed impuri- 
ties (corresponding to 5% density). By considering all 
possible abstract graph cmbcddings in the finite system, 
we calculate exact 1/T series for the local susceptibilities 
in the sample. Analyzing these as before, we obtain re- 
liable quantitative results down to about T = 0.5 J. The 
local susceptibility histograms for several temperatures 
are shown in Figure [D We also calculate the suscepti- 
bilities at the oxygen atom sites that the experiments 
measure; the corresponding histograms are displayed in 
Fig.H 

The sample is sufficiently large that the obtained his- 
tograms are accurate representations of an infinite sys- 
tem; in particular, all features seen in Fig.[2]are real. The 
peaks are not delta functions because all sites feel the im- 
purities, and the different peaks correspond to different 
locations relative to one or several impurities. We esti- 
mated the errors by dividing the sample into ten pieces, 
getting the histogram for each, and calculating deviations 
from the average. The resulting relative error on the y 
axis in Fig. [2] is roughly 5%. There is also a systematic 
error along the x axis coming from the truncated series 
and the Fade approximants, but from the single impurity 
study we believe the results arc reliable for the tempera- 
tures shown. 

To the first approximation, we can say that the are two 
major groups of sites: one is the sites next to impurities 
(susceptibilities in the histogram near xi), and the other 
is the rest of the sites (susceptibilities near Xbuik)- In the 
second group, we can also clearly see a subgroup formed 
by the sites that are second neighbors to some impurity 
and have susceptibilities near X2) X3- 

Turning to the ^^O NMR experiments, each oxygen can 
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FIG. 2: Histogram of local susceptibilities at the sites of 
Kagome lattice in the presence of 5% of nonmagnetic im- 
purities. Fade approximant [6,6] was used to analyze exact 
series at each site of a large but finite 60 x 60 x 3 sample with 
randomly placed impurities, x is in units (g/iB)^ /J- 



be associated with a bond < ij > of the Kagome lattice 
and feels only these two sites i and j. The appropriate 
histogram of such Xi + Xj is in Figure [31 

In the first approximation, we sec four major groups of 
peaks which can be roughly identified as follows. Imag- 
ine that only the susceptibility of the sites next to the 
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FIG. 3: Histogram of the susceptibilities at oxygen sites in 
the presence of 5% of nonmagnetic impurities. The system is 
the same as in Fig.[2]and we used the same Fade approximant 
[6,6]. 



impurities is modified and is xi j while the rest of the Cu 
sites have Xbuik (in particular we join the X2 and X3 with 
the Xbuik, even though they are still separately visible in 
Fig. [2). The leftmost peak in Fig. [3] is located near xi 
and comes from the oxygens next to impurity that feel 
only one Cu site, and so the value is roughly half com- 
pared to all other oxygens. The dominant peak is near 
2Xbuik and comes from oxygens that are in the bulk - nei- 
ther near impurity, nor near Cu next to impurity. The 
rightmost peak is located around 2xi and comes from 
oxygens that are interacting with two Cu sites that are 
next to impurity. Finally, the group Xi -|- Xbuik comes 
from oxygens that are interacting with one site near and 
one not near impurity. Again, the finer features and the 
finite widths of the peaks originate from the variation of 
the susceptibility with the distance from various impuri- 
ties a given site feels. The volume fractions of the four 
groups of peaks in Fig. [3] say at T = 0.6J, are, from 
the left, roughly 8%, 70%, 15%, 5%. If we calculate how 
many oxygens of each type other than 2xbuik are near 
a single impurity and then multiply it by the density of 
impurities and divide by the ratio of the number of oxy- 
gens to the number of Kagome sites (two), we would get 
the following percentages 10%, 65%, 20%, 5%; these are 
slightly different from the exact numbers because some 
sites are close to more than one impurity. 

Thinking about the experiments, the NMR lines are 
too broad to resolve the groups 2xbuik, Xi + Xbuik, and 
2xi, but just enough to separate these from the group 



Xi- We hope that our study will stimulate further NMR 
experiments with more site resolution. 



III. DIRAC SPIN LIQUID: MEAN FIELD WITH 
IMPURITIES 

In this section we explore the Kagome system with 
non-magnetic impurities at low temperatures within 
a particular spin liquid proposal of Hastings^i and 
Ran et al^. In the clean system, one has fermionic 
spinous hopping on the Kagome lattice in the presence 
of TT flux through every hexagon. Ran et al^ showed 
that the corresponding spin wavefunction obtained by 
Gutzwiller projection has very good variational energy.— 
Again we consider the cases of a single impurity and 5% 
density of impurities. We take a crude model of the spin 
liquid in the presence of impurities by simply removing 
the vacancy sites (and the links) in the spinon hopping 
problem. We will argue that while this is not fully ade- 
quate, particularly for detailed features near impurities, 
the results for the overall behavior of the bulk of the lo- 
cal susceptibility histograms are robust. We will further 
discuss possible improvements of our approximation in 
Sec. Himi 

In the mean field the fermions are free. The spec- 
trum of the clean system has a flat band at the low- 
est energy containing 1/3 of states with gap to the next 
band. The system has Dirac points at fillings 1/2, 2/3 
and 5/6. In our approach, adding impurities does not 
destroy the flat band (which can be expressed in terms 
of localized states), but also produces localized states in 
the gap. For a flnite concentration of impurities, which 
is our main interest, the Dirac points are filled with fi- 
nite density of states (DOS). We consider the system at 
half filling which corresponds to spin 1/2. The pure sys- 
tem has Dirac nodes and a linearly vanishing density of 
states at the Fermi level, while the impurities induce a 
finite density of states, and one of the main questions we 
address is how it is distributed locally. 

If {V'n(*)} is the set of single-particle wavcfunctions, it 
is easy to show that the local susceptibility at tempera- 
ture T is given by 

Xloc(^) = ^^El'/'"«lV(en)(l-/(en)) (2) 

n 

where /(e) ~ l/(e(^"'^'/-^ + 1) is the Fermi distribution. 
For a given configuration of impurities, we obtain the 
single-particle orbitals by numerical diagonalization and 
use this formula to calculate the local susceptibility. In 
the T limit, the local susceptibility is simply pro- 
portional to the local density of states at the Fermi level, 

t'loc(«) = E„ \tpn{i)\'^S{£n ~ M)- 

Let us briefly describe the localization properties of 
the single-particle states in our system with 5% density 
of impurities. We expect that all states are eventually 
localized, since this is a time-reversal invariant Anderson 
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localization problem in two dimensions with no special 
symmetries. The localization is weaker in the spectral 
regions where the density of states is high and is stronger 
where the DOS is low, e.g., in the vicinity of the Dirac 
nodes of the pure spectrum. However, even for the latter 
as happens around half-filling, the localization lengths 
are still very large. For example, wc calculate the partic- 
ipation ratios of the states near half-filling and find that 
the corresponding volume fraction of participating sites 
decreases only from about 14% in a 30 x 30 x 3 system 
to about 10% in a 60 x 60 x 3 system, indicating that 
these states are still extended over our system sizes up to 
60 X 60 x 3. On the other hand, such long distance local- 
ization physics does not affect our measurements of local 
properties at finite temperatures presented here. For ex- 
ample, even a single disorder sample in a system as small 
as 10 X 10 X 3 produces the susceptibility histogram which 
is essentially the same as from much larger samples. 

Some explanations are in order before we show the re- 
sults for the local susceptibilities. Throughout, we keep 
the spinon hopping t fixed and vary the temperature. 
All presented susceptibilities are in units of {q^.b)" /t- In 
a more realistic calculation, the spinon hopping would 
need to be found self-consistently for a given Heisenberg 
exchange J and temperature T. In such treatments of 
pure systems, one typically finds that the self-consistent 
t vanishes above some temperature of order J (e.g., in 
the rcnormalizcd mean field scheme this temperature is 
0.75 J). There is no sense of speaking about spinous 
above this temperature, and the system is more appropri- 
ately described in terms of weakly correlated individual 
spins. On the other hand, t becomes nonzero at lower 
temperatures, suggesting that the system enters the cor- 
related paramagnet regime. Below the onset tempera- 
ture, the self-consistent t quickly reaches the values simi- 
lar to the zero-temperature limit. It is this regime where 
t{T) ~ t{T = 0) that we are describing. 

We can estimate the spinon hopping amplitude from 
self-consistent mean field conditions. For example, in the 
so-called renormalized mean field scheme, this reads 



A. Pure System and System with a Single Impurity 

As mentioned, the local susceptibility near an impu- 
rity depends strongly on how the spin liquid is modified 
near this impurity. However, for illustration, we show 
the local susceptibility of the first few neighbors for the 
case where we simply remove the appropriate links from 
the fermion hopping problem. For the size 44 x 44 x 3 
this is shown in Fig. U) The susceptibilities vanish lin- 
early when T —t Q, since the pure system is perturbed 
only by a single impurity here. In our approximation, 
the susceptibility xi of the site next to the impurity is 
larger than the bulk susceptibility. This happens because 
we take all bonds to have the same strength, while the 
Xi sites have fewer bonds. This is an example where 
we believe the simplified treatment near the impurity is 
inadequate and the xi results should not be taken seri- 
ously. A more accurate self-consistent treatment would 
likely require a stronger hopping amplitude connecting 
two such sites next to an impurity, and this would de- 
crease xi- We discuss such possible improvements over 
the current treatment in Sec. IIII Dl 
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FIG. 4: [Color Online] The susceptibilities of first few neigh- 
bors near single impurity in the spin liquid with tt fluxes 
through each hexagon on Kagome lattice, x is in units 

{9t^B?/t. 



where Qs is some renormalization factor (which can be 
found, e.g., by comparing measurements before and after 
the projection). In the case with tt fiux through each 
hexagon and at T = 0, Hastings found K/jo-Ziir)! = 
2 • 0.221. Using a reasonable gs = 4, we then find 
t = 0.66J. On the other hand, from the study of the 
projected wavefunctions. Ran et al^ suggest 

t « 0.4J « 70K. (4) 

These are crude but reasonable estimates of the spinon 
band energy scales. 



We also plot the shape of local susceptibility near a 
single impurity in Figure [5] The "ringing" that we see 
is a more universal aspect of the Knight shift texture 
around the impurity. The corresponding wavevectors are 
the mid-points of the Brilloin zone edges and arc charac- 
teristic for this spin liquid. One can obtain some analytic 
understanding of the mean field results by working per- 
turbatively in the non-magnetic impurity strength. One 
finds that (xioc — x)/x decays away from the impurity 
with a 1/r envelope on length scales 1 ^ r <C v/T and 
exponentially on larger scales (here v is the Fermi veloc- 
ity, and wc referenced xioc by x remembering that the 
susceptibilities vanish in the T — s- limit). 
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FIG. 5: The susceptibilities near single impurity in the spin 
liquid with tt fluxes through each hexagon on Kagome lattice 
at temperature T = 0.054f where t is the hopping amplitude. 
In this figure, the color of a given point in the plane is given 
by the value at the nearest Kagome site. 



B. 5% Density of Impurities 

For the system with 5% density of impurities, the his- 
togram of local susceptibilities at Kagome (Cu) sites is 
in Figure [5] and at oxygen sites in Figure [71 
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FIG. 6: The histogram of susceptibilities on the Kagome lat- 
tice sites in the spin liquid of fermionic spinous with n fluxes 
through hexagons in the presence of 5% of impurities, treated 
by simply removing the sites from the spinon hopping prob- 
lem. The histogram is the average of 300 samples of size 
20 X 20 X 3. X is in units {gfj.B)^/t. 

We can understand the main features in the xioc his- 
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FIG. 7: The histogram of susceptibilities on oxygen sites in 
the spin liquid in the presence of 5% of impurities. The sys- 
tems is the same as in Fig. [S] The histogram is the average 
of 300 samples of size 20 x 20 x 3. 



tograms as follows. At the highest temperature shown, 
T = 0.44t, there are two main groups, just as in the high 
temperature series analysis. The group with the higher 
X comes from the sites that are next to the impurities, 
while the other with the lower x comes from the rest 
of the sites (bulk) and contains the main weight of the 
sites. As the temperature is lowered, the two features 
move towards smaller susceptibilities, but the bulk one 
moves faster. Also, between T = OAt and T = 0.2t, 
we see a new feature splitting away from the bulk and 
moving slower; we can roughly identify this feature with 
the sites that have impurities as their second neighbors. 
The continuing marching of the bulk sites that are fur- 
ther from the impurities to smaller x is what the clean 
Kagome lattice would do because of the vanishing den- 
sity of states, see Fig. |4l As we lower the temperature, 
however, the influence of the impurities spreads to fur- 
ther and further neighbors, and the local susceptibility 
at the sites that already feel the impurities starts to slow 
down and eventually saturates at a finite local x- 

At T = O.lt, we already see no sharp features associ- 
ated with specific locations away from an impurity since 
all the sites already feel many impurities. The broad 
peak corresponds to the "bulk" sites that are furthest 
from impurities; it remains visible since this group of 
sites contains the largest volume fraction among the suc- 
cessive "layers" around impurities. The bulk peak stops 
its motion below T = O.lt; its eventual location cor- 
responds roughly to the typical impurity-induced local 
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density of states at these bulk sites. Note that this value 
Xpeak ~ 0.02/t is about two times smaller than the aver- 
age over all sites x = 0.0425/t; the latter contains con- 
tributions from the sites close to impurities where the 
induced density of states is typically larger. 

This behavior of the bulk peak is our main and robust 
prediction. On the other hand, the behavior of the other 
peaks visible at high temperatures depends on specific 
details of the treatment of the impurity; these features 
should not be trusted in our simplified approach and we 
will discuss the ways to improve the treatment in Sec- 
tion [mDl 

Thinking about the NMR experiments, the bulk peak 
susceptibility at low temperatures Xpoak ~ 
0.04/J, where we used Eq. This is about one half 
of the bulk spin susceptibility at T = 2 J from the high 
temperature series study, Fig.[T] In the ^^O NMR exper- 
iments, this ratio is between one half and one third, so 
our spin liquid results are reasonable. 

To conclude this section, let us see if the disorder- 
induced density of states can give sensible thermody- 
namic properties in the Herbertsmithite. In our treat- 
ment with 5% impurity concentration, the DOS at the 
Fermi level is 



i^iep) = 0.085/t 



(5) 



per site per single spin species. This would give average 
spin susceptibility of the Cu in the Kagome layer 



X = 2K6f)mI = 6.4.10-^x 



mol Cu 



(6) 



One should not take these numbers literally, since i^iep) 
depends quantitatively on the treatment near impurities, 
which is very crude here. Still, an estimate of t like that 
in Eq. (g]) produces x ~ 10"'^ emu/mol Cu that compares 
reasonably with the experimental estimate of the intrinsic 
Kagomc layer susceptibility in Ref. [^. Furthermore, in 
the mean field treatment, spinous will contribute to the 
specific heat as 

^ , N , 9 ^ . „ fc^r Joule 
e.=2.M-4T = 4.6x^ ' 

which again yields reasonable numbers for temperatures 
of several Kelviniii^ 



C. Gutzwiller Wavefunction Study of the Local 
Magnetizations 

In this section, we pursue the spin liquid proposal 
one step beyond the mean field. For the slave fermion 
approaches, one can achieve this by performing the 
Gutzwiller projection of the singlet mean field state into 
the space with precisely one fermion per site. This ap- 
proach is appealing since one obtains a physical spin 
wavefunction. which can be viewed as a variational state 



for an energetics study (and one can also study proper- 
ties of this spin state). As we have already mentioned. 
Ran et al^ performed the energetics study of the clean 
Kagome system and found the discussed 7r-hexagon spin 
liquid to have the lowest trial energy. We expect that the 
same construction on a lattice with vacancies will have 
good energetics also for small concentration of impurities. 
It is interesting to explore this issue in more detail in the 
future, and we comment on this in Sec. lIII Dl below, while 
here we want to take a crude look whether such proposal 
can explain the observed inhomogeneous Knight shifts in 
the Herbertsmithite^ii The specific results that we find 
here are as follows. The spin susceptibilities are always 
positive in the mean field. We show that they can become 
negative and the Knight shift distributions are actually 
broader when we go beyond the mean field. 

Unfortunately, one can not access the local suscepti- 
bilities from the study of the singlet ground state. To 
circumvent this, we consider trial states with non-zero 
^t^oj and study how the total spin is distributed among 
the lattice sites. The idea is that the pattern of the lo- 
cal magnetizations m| = (S'f ) in such weakly polarized 
states is similar to that of the local susceptibilities at 
low temperatures. This approach has been used in ex- 
act diagonalization studies of spin systems with impuri- 
ties, where one considers mf in the lowest-energy states 
in the sectors with nonzero S^^^. We are trying this in 
the variational spin liquid studies in the present Kagome 
problem and also for the organic spin liquid on the tri- 
angular latticci^ (We have also tried this approach in a 
one-dimensional chain with vacancies^ and reproduced 
the non-trivial qualitative behavior predicted by Eggert 
and Afiieck.l&) 

Consider a mean field state with different chemical po- 
tentials for the two spin species, /X| < fxp The orbitals 
with /ij^ < e < ^1 are occupied by the t spins but not by 
the I spins, which produces non-zero S^^^ = {N-^—Ni)/2. 
In the mean field, the local magnetization is given by 



(8) 



Observe the similarity of this expression to that for the lo- 
cal susceptibility at finite temperature, Eq. The dif- 
ference is that here the contributing orbitals are selected 
by a sharp energy window while in Eq. ([2|) the 

window is soft proportional to (—df/de) with width set 
by the temperature. If there is a finite density of states 
at the Fermi level and if the properties of the orbitals 
do not change strongly with the energy, the two win- 
dows of similar width should give similar results for the 
local magnetization or susceptibility normalized by the 
corresponding average values. The above two windows 
roughly match when /Lt| — /i| « 4T. In our numerical 
calculations, we indeed observe such correspondence of 
the mioc and xioc distributions. 

From the mean field state with iioii-zcro S^Q^^ WG con- 
struct the spin wave function by Gutzwiller projection. 
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Kagome 10x10x3 with 15 vacancies; N^=146, N^=139 Kagome 10x10x3 with 15 vacancies; N^=146, N^=139 




FIG. 8; [Color Online] Distribution of the local magnetization 
in the vr-flux spin liquid state with non-zero ^tot • The sample 
is 10 X 10 X 3 with 15 randomly placed vacancies; the total 
magnetization is Sfot = (-^T ~ -^i)/2 ~ 7/2. We compare the 
mean field and the Gutzwiller-projected results for exactly 
the same system and find that the projection broadens the 
histogram by about a factor of 1.5-2, which can be clearly 
seen despite the noisiness of the data in the single random 
sample. 



Note that S^^^ is preserved by the projection, so the aver- 
age magnetization is unchanged. However, the local mag- 
netizations are modified. Our observations in the present 
system can be summarized by the folloviring crude rule: 
The deviation of the site magnetization from the average 
is enhanced by the projection by a typical factor between 
1.5 and 2. This implies that the distribution of the lo- 
cal magnetization is broadened by the same factor upon 
the projection. Figure [8] illustrates this for a 10 x 10 x 3 
sample with 15 vacancies and — Ni = 7. 

We have tried several lattice sizes, different disorder 
samples, and different magnetizations S^^^ = (iV| — 
Ni)/2, and in all cases observed the above rule. This 
rule is not strict but holds for a majority of sites. 
More precisely, the histogram of the ratios {m^^^™ — 
"rh) / {m™^^'^^'^^^ — ifi) is peaked near 1.5 — 2, which is the 
claimed typical enhancement factor for the deviations 
from the average. In particular, while the mean field 
local magnetization is always non-negative, sec Eq. ([5]), 
the local magnetization can become negative after the 
projection. It can be opposite to the bulk magnetization 
because of the antiferromagnetic correlations present in 
the system. 

Let us say few words about our choice of , Ni in 
Fig. [HI On one hand, if we take — = 1, we 
worry if the results are dominated by the specific single- 
particle state that happened to lie at the Fermi level. 
This is a concern because of our view that all states 
will be eventually Anderson-localized. However, by con- 
sidering the participation ratio for the states near the 
Fermi level in the present sample, we find that they are 
still spread over about one third of the sites. We have 



FIG. 9: [Color Online] Oxygen NMR local fields correspond- 
ing to the same system as in Fig. |8l 



checked that even for such small 7V| — the results 
for the broadening of the distribution by the projection 
arc qualitatively similar. On the other hand, we do not 
want N-f — Ni to be large, which could take us far from 
the region in the spectrum near the Fermi level and also 
would correspond to physically inaccessible magnetiza- 
tion. Our N-^ = 146, Ni = 139 is a compromise, since 
a few typical states around the Fermi level are sampled, 
and also since it corresponds to about 2.5% of the max- 
imal polarization, which is reasonable also from the ex- 
perimental perspective. For example, using the estimate^ 
Xintrinsic = 1.25 • 10"'^ cmu/mol Cu, the field of 7 Tesla 
in the ^^O experiments^ magnetizes the system to about 
1.5% of the maximal polarization of the Kagome layer. 

To conclude this section, we find that the Gutzwiller 
trial wavefunction treatment docs not modify the mean 
field results qualitatively but widens the magnetization 
distributions by about a factor of twoi^ In particular, at 
low temperature the distribution will have a sizable wing 
on the negative values, which is consistent with the NMR 
experiments!^ To describe the ^^O NMR line, we need to 
histogram the field rrii + mj associated with bonds < ij > 
of the Kagome lattice. This is shown in Fig.[9l The wing 
at negative fields is weaker than in the Cu NMR because 
the oxygens perform some "averaging" of the local 
(while the total magnetization is positive). We reiterate 
here the crudeness of our study; in particular, we are not 
attempting here to explain the actual experimental his- 
togram in the ^^O NMR, which looks nearly symmetric 
with respect to the peak at the lowest temperature. Our 
main message is that within the spin liquid approach, by 
going beyond the mean field, we expect broader distribu- 
tion of the Knight shifts with both positive and negative 
susceptibilities due to the underlying correlations. 
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D. Issues with the Spin Liquid Treatment and 
Possible Improvements 

The main reservation about the presented spin hquid 
study is that we modeled the effect of impurities by sim- 
ply removing the Zn sites in the spinon mean field. We 
did not perform a self-consistent mean field or an ener- 
getics study of the system with impurities. However, we 
did think how such more accurate treatments could affect 
our results and tested the robustness of the features by 
trying different modifications that such treatments might 
require. 

One minor concern is that even though the system 
is inhomogeneous, we require the half- filling for spinous 
only on average over the whole sample and not at each 
site. The latter is likely a better treatment since then the 
Gutzwiller projection is less drastic and the mean field is 
closer to the projected state. Wc find that in our samples 
the deviations of the local densities from the average are 
small and smaller in relative terms than the deviations 
of the local susceptibilities that wc arc interested in. We 
also tested small systems achieving the half-filling pre- 
cisely at all sites by adjusting local chemical potential 
and found that the susceptibilities did not change qual- 
itatively or significantly quantitatively and only at the 
the lowest temperatures the histograms became some- 
what broader. 

A more serious concern is that wc took the spinon hop- 
ping amplitudes to be the same as in the clean system, 
instead of finding them self-consistently. Thus, we ex- 
pect that at low temperatures the hopping connecting 
two sites that are both next to an impurity should be 
larger than for other bonds. Consider, for example, the 
renormalized mean field scheme, Eq. ([3]). When we calcu- 
late the bond expectation value Qj = ifjfi) in our mean 
field starting with all Uj equal, which can be viewed as 
a first step in the self-consistency iteration, we find that 
it is the largest on the bonds next to impurities. We 
have not pursued this scheme any further since in the 
present case it would actually lead to a dimer state in 
the clean system, which we know is a pathology of such 
treatment when applied to the spin-1/2 system. Roughly, 
one should not be using the same renormalization factor 
gs irrespective of the bond expectation value. By com- 
paring with the Gutzwillcr-projected measurements, one 
typically uses g,, = 4 for translationally invariant states, 
but one should take a smaller Qs = 2 for the dimer states. 
One could use a different mean field scheme to capture 
this, e.g., the scheme due to Hsuf^i^ which leads to the 
self-consistency conditions of the form 

, _ 3J,,C,(l-16|G,n .q. 
(1 + 1610,1^)2 ■ 

This suppresses the dimerization tendency of the less 
elaborate mean field but may be doing it already too 
much, while we may want to retain some possibility of 
local dimcrs in the system with impurities. We leave a 
detailed exploration of such approaches and whether they 



can capture the actual energetics in the system with im- 
purities to future work. 

It is interesting to notice however that if wc took the 
most naive treatment Fig. 31 at low enough temperatures 
the nearest neighbor susceptibility is larger than twice 
the bulk and a little peak starts to emerge on the right 
of the bulk peak in Fig. [7] for T « 0.1 It. If this were 
the case then the peaks in Ref. at low temperatures 
would have the opposite interpretation - bulk peak would 
become closer to zero and the next to impurity oxygen 
peak further from zero. However we cannot make this 
conclusion at this point since the bonding of the spins 
next to the impurity will decrease the local susceptibility 
near the impurity. We only would like to point out that 
one should be careful when interpreting the experiment. 

For the treatment of the spin liquid near impurity, we 
only note some crude things that we tried to see how 
our results might be affected. Motivated by the above 
observation, we have studied local susceptibilities at low 
temperatures when we multiplied by 2 all the bonds con- 
necting the two sites that are both next to an impurity. 
Certainly, this has a strong effect on the susceptibility xi 
of these sites - it in fact becomes smaller than the bulk 
one, which is opposite to the results in the preceding sec- 
tion where the xi is the largest. Thus, such treatment 
can make the results for xi niore in line with the series 
study. Sec. [Tll and exact diagonalization study, Ref. [33l . 
that suggest that this susceptibility is depressed at low 
temperatures. On the other hand, the above change of 
the hopping strengths near impurities does not have sig- 
nificant effect on the bulk susceptibility histogram peak 
and its marching to lower values with decreasing temper- 
ature (except, of course, for small numerical differences). 
We have also tried modifying not only the strengths of 
the bonds but also fluxes next to impurities. For ex- 
ample, if we view the proposed 7r-hexagon spin liquid as 
a time reversal invariant way of performing flux attach- 
ment and flux-smearing mean field, obtaining Dirac spin 
liquid instead of a chiral spin liquid (see discussion in 
Sec. lie of Ref. |4^ . it suggests that we should also re- 
move a TT flux when we remove a spin. In an exploratory 
Gutzwiller energetics study with impurities, we indeed 
find that this often improves the trial energy. As far as 
we are concerned here with the local susceptibilities, such 
local flux modification of course has an effect on the sites 
near the impurities, but has little effect on the described 
qualitative behavior of the bulk peak. 

We also remark that making all hexagon fluxes equal to 
zero^l results in a qualitatively different behavior of the 
bulk of the susceptibility histogram. Here the clean sys- 
tem has a finite density of states at the Fermi level and 
therefore a nonzero susceptibility at low temperatures. 
We tried this for the system with 5% density of impuri- 
ties and indeed found that the bulk peak doesn't march 
down. Instead, it approaches a relatively large value of 
susceptibility and broadens to both high and low values 
of susceptibility. We then propose that the experiments 
suggest more a Dirac-like spin liquid than the one with 
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a full Fermi surface. 



IV. CONCLUSIONS 

Using the high temperature series expansion, we calcu- 
lated local susceptibility for the nearest neighbor Heisen- 
berg antiferromagnet on Kagome lattice down to about 
T w 0.4J. The resulting histogram of susceptibili- 
ties can be directly compared to the experiment to test 
whether the Heisenberg model with impurities applies to 
ZnCu3(OH)6Cl2. The calculated histogram contains a 
lot of fine structure some of which we hope can be ob- 
served if the resolution of the experiments improves. One 
issue not considered in experiments is that the local sus- 
ceptibility xi next to impurities is larger than the bulk 
value by as much as 20% at high temperatures between 
0.5 J and 2 J. This should affect the peak associated with 
the oxygens close to impurity, and one should be care- 
ful at interpreting the observed features at the interme- 
diate and lower temperatures. Our work suggests that 
Xi drops below the bulk value at low temperatures and 
is consistent with earlier exact diagonalization study^ 
which found that these sites tend to form singlets. The 
latter was invoked in Ref. to explain a sharp feature in 
the NMR signal appearing near low temperatures. This 
is appealing, but more work needs to be done to under- 
stand the origin of this feature, particularly in the pres- 
ence of the sizable impurity concentration, and also with 
the added bulk peak contribution that we expect in the 
spin liquid regime. 

At low temperatures we assumed that the system forms 
a spin liquid of Ran et al^ We calculated the susceptibil- 
ity histograms in the mean field theory and went beyond 
mean field by using Gutzwiller projection. Our results 
for the local susceptibility should not be trusted quan- 
titatively close to impurity, since these depend on the 
details of the treatment of the spin liquid near the impu- 
rities, while we took a crude approach and only outlined 
how one can make it more accurate. In particular, we 
are not able to discuss histogram features associated with 
sites close to impurities. However we have a robust pre- 
diction for the behavior of the peak coming from the bulk 
of the sites (Fig. [6]) : The bulk peak marches down with 
decreasing temperature; its weight is gradually spread as 
further layers of sites around the impurities start to feel 
their share of the impurity-induced density of states; the 



broadening saturates when all sites start to feel several 
and more impurities. We expect this qualitative descrip- 
tion to be insensitive to the details of the mean field in 
the immediate vicinity of the impurities. It is sufficient 
to know that the clean spin liquid has Dirac density of 
states and that the bulk sites do not see the impurities 
down to some temperature scale, while any treatment 
near the impurities represents some local disorder that 
always induces some density of states at the bulk sites. 
The effect of the Gutzwiller projection is crudely to in- 
crease the variation of the local susceptibility by about a 
factor of two. In particular, as the peak spreads some of 
the susceptibilities become negative. 

This behavior is consistent with the observed Knight 
shifts in the ZnCu3(OH)6Cl2 experiments. Such behav- 
ior observed in the presence of a sizable density of im- 
purities then suggests that the clean spin liquid would 
have a vanishing density of states at the Fermi level. 
Thus it distinguishes this spin liquid from the proposal 
with a spinon Fermi surface for which the peak would 
not march down but instead saturate at a large positive 
value. On the other hand, since the Algebraic Vortex 
Liquid of Ryu et al^ has Dirac-like gapless character, it 
would likely produce inhomogeneous Knight shift behav- 
ior similar to the Dirac spin liquid. 

In the future work it would be good to treat the spin 
liquid right near impurities more accurately, e.g., in a 
self-consistent mean field or in an energetics study, and 
consequently make a reliable prediction for the spin sus- 
ceptibility at these sites. Quantitative comparison with 
other approaches such as exact diagonalization would 
also be useful. It would be interesting to consider quan- 
titatively the gauge theory description of such spin liq- 
uid with disorder-induced finite density of states, though 
we do not expect our mean field results to be modified 
significantly. In a companion work,— we are studying 
non-magnetic impurities in the triangular lattice antifer- 
romagnet and explore spin liquid with a spinon Fermi 
surface, which may be relevant for the organic spin liq- 
uid material k-(ET)2Cu2(CN)3. 
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